Acute Changes in CEWL as a Function of Temperature: Data Analysis

Calvin Davis, Savannah Weaver

2022

Packages

if (!require("tidyverse")) install.packages("tidyverse") 
library("tidyverse") #tidyr, ggplot, dplyr, etc
if (!require("ggpubr")) install.packages("ggpubr") 
library("ggpubr")
if (!require("lme4")) install.packages("lme4") 
library("lme4") # mixed-effect models
if (!require("lmerTest")) install.packages("lmerTest") 
library("lmerTest") # p-values
if (!require("car")) install.packages("car") 
library("car") # test for VIF
if (!require("broom.mixed")) install.packages("broom.mixed")
library("broom.mixed") # model export
if (!require("emmeans")) install.packages('emmeans')
library('emmeans')
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # model export
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown format

Background

This code walks through the statistical models and figures for a study on the rate of change of cutaneous evaporative water loss (CEWL) in response to temperature change in Western Fence Lizards (Sceloporus occidentalis) (published in the Journal __ in 2023). Data was collected and analyzed by Calvin Davis and Savannah Weaver, under the advising of Dr. Emily Taylor at California Polytechnic State University.

Load Data

Get the RDS files created in the data wrangling Rmd:

temp_time_series <- read_rds("Data/temp_time_series.RDS")
CEWL_time_series <- read_rds("Data/CEWL_time_series.RDS")

Stats & Models

Temp Correlation

Check the correlation between dorsal and cloacal temps:

temp_corr <- lm(data = temp_time_series, 
                calibrated_dorsal_temp ~ calibrated_cloacal_temp)
summary(temp_corr)
## 
## Call:
## lm(formula = calibrated_dorsal_temp ~ calibrated_cloacal_temp, 
##     data = temp_time_series)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0389 -0.7784  0.0075  0.5911  5.3351 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             6.966310   0.049820   139.8   <2e-16 ***
## calibrated_cloacal_temp 0.747869   0.001899   393.8   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.388 on 25249 degrees of freedom
## Multiple R-squared:   0.86,  Adjusted R-squared:  0.8599 
## F-statistic: 1.55e+05 on 1 and 25249 DF,  p-value: < 2.2e-16

They are actually pretty different, so run parallel models of CEWL, one with dorsal temp, and one with cloacal temp.

CEWL (raw) ~ temp LMM

First, test linear effects of temperature:

CEWL_LMM_02a <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             calibrated_cloacal_temp *
                             treatment +
                             (1|date/lizard_ID))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues
CEWL_LMM_02b <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             calibrated_dorsal_temp *
                             treatment +
                             (1|date/lizard_ID))

Compare to the inclusion of polynomial factors:

CEWL_LMM_03a <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             poly(calibrated_cloacal_temp, 2)*treatment +
                             (1|date/lizard_ID))

anova(CEWL_LMM_03a, CEWL_LMM_02a)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_02a: CEWL_g_m2h ~ calibrated_cloacal_temp * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_02a    9 87929 88002 -43956    87911                         
## CEWL_LMM_03a   12 81676 81774 -40826    81652 6259.3  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
CEWL_LMM_03b <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             poly(calibrated_dorsal_temp, 2)*treatment +
                             (1|date/lizard_ID))

anova(CEWL_LMM_03b, CEWL_LMM_02b)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_02b: CEWL_g_m2h ~ calibrated_dorsal_temp * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_02b    9 88775 88848 -44379    88757                         
## CEWL_LMM_03b   12 83092 83190 -41534    83068 5689.1  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The poly models have significantly better explanatory power than the linear models.

Also try log transformation to make sure a polynomial effect is best:

CEWL_LMM_04a <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             (log(calibrated_cloacal_temp))*treatment +
                             (1|date/lizard_ID))

anova(CEWL_LMM_04a, CEWL_LMM_02a)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_04a: CEWL_g_m2h ~ (log(calibrated_cloacal_temp)) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_02a: CEWL_g_m2h ~ calibrated_cloacal_temp * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)
## CEWL_LMM_04a    9 88930 89003 -44456    88912                     
## CEWL_LMM_02a    9 87929 88002 -43956    87911 1000.4  0
anova(CEWL_LMM_04a, CEWL_LMM_03a)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_04a: CEWL_g_m2h ~ (log(calibrated_cloacal_temp)) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_04a    9 88930 89003 -44456    88912                         
## CEWL_LMM_03a   12 81676 81774 -40826    81652 7259.7  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
CEWL_LMM_04b <- lme4::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             log(calibrated_dorsal_temp)*treatment +
                             treatment +
                             (1|date/lizard_ID))

anova(CEWL_LMM_04b, CEWL_LMM_02b)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_04b: CEWL_g_m2h ~ log(calibrated_dorsal_temp) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_02b: CEWL_g_m2h ~ calibrated_dorsal_temp * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)
## CEWL_LMM_04b    9 89762 89835 -44872    89744                     
## CEWL_LMM_02b    9 88775 88848 -44379    88757 987.18  0
anova(CEWL_LMM_04b, CEWL_LMM_03b)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## CEWL_LMM_04b: CEWL_g_m2h ~ log(calibrated_dorsal_temp) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## CEWL_LMM_04b    9 89762 89835 -44872    89744                         
## CEWL_LMM_03b   12 83092 83190 -41534    83068 6676.3  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The logarithmic models are NOT better than the linear or polynomial ones.

Check assumptions:

plot(CEWL_LMM_03a)

plot(CEWL_LMM_03b)

vif(CEWL_LMM_03a)
##                                                    GVIF Df GVIF^(1/(2*Df))
## poly(calibrated_cloacal_temp, 2)           9.312576e+05  2       31.064721
## treatment                                  1.100911e+00  2        1.024326
## poly(calibrated_cloacal_temp, 2):treatment 9.312225e+05  4        5.573547
vif(CEWL_LMM_03b)
##                                                   GVIF Df GVIF^(1/(2*Df))
## poly(calibrated_dorsal_temp, 2)           5.528187e+05  2       27.267523
## treatment                                 1.107683e+00  2        1.025897
## poly(calibrated_dorsal_temp, 2):treatment 5.527949e+05  4        5.221803

They both have a really weird spike, but otherwise there’s no fanning or different pattern. VIF is negligible.

Re-run with p-values:

CEWL_LMM_besta <- lmerTest::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             poly(calibrated_cloacal_temp, 2) * treatment +
                             (1|date/lizard_ID))
summary(CEWL_LMM_besta)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(calibrated_cloacal_temp, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 81611.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6892 -0.6401 -0.0881  0.5378  7.8576 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 11.773   3.431   
##  date           (Intercept) 23.794   4.878   
##  Residual                    1.468   1.212   
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                                      Estimate Std. Error
## (Intercept)                                            0.6557     2.4503
## poly(calibrated_cloacal_temp, 2)1                   1422.2807    62.9113
## poly(calibrated_cloacal_temp, 2)2                  -3419.1801    82.4963
## treatmentCooling                                       6.6956     1.5010
## treatmentHeating                                      16.7664     1.5897
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling -1162.0238    62.9471
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling  3565.6587    82.5285
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating -1223.6577    62.9510
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating  3338.2055    82.5256
##                                                            df t value Pr(>|t|)
## (Intercept)                                            5.4550   0.268 0.798840
## poly(calibrated_cloacal_temp, 2)1                  25235.9702  22.608  < 2e-16
## poly(calibrated_cloacal_temp, 2)2                  25238.1540 -41.446  < 2e-16
## treatmentCooling                                      30.3218   4.461 0.000104
## treatmentHeating                                      29.8965  10.547 1.36e-11
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling 25235.9622 -18.460  < 2e-16
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling 25238.1521  43.205  < 2e-16
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating 25235.9530 -19.438  < 2e-16
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating 25238.1526  40.451  < 2e-16
##                                                       
## (Intercept)                                           
## poly(calibrated_cloacal_temp, 2)1                  ***
## poly(calibrated_cloacal_temp, 2)2                  ***
## treatmentCooling                                   ***
## treatmentHeating                                   ***
## poly(calibrated_cloacal_temp, 2)1:treatmentCooling ***
## poly(calibrated_cloacal_temp, 2)2:treatmentCooling ***
## poly(calibrated_cloacal_temp, 2)1:treatmentHeating ***
## poly(calibrated_cloacal_temp, 2)2:treatmentHeating ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.159                                                        
## ply(c__,2)2  0.164 -0.914                                                 
## tretmntClng -0.336  0.260    -0.268                                       
## tretmntHtng -0.324  0.244    -0.251     0.505                             
## pl(__,2)1:C  0.159 -0.999     0.914    -0.259 -0.244                      
## pl(__,2)2:C -0.164  0.914    -1.000     0.268  0.251 -0.913               
## pl(__,2)1:H  0.159 -0.999     0.914    -0.259 -0.244  0.999     -0.913    
## pl(__,2)2:H -0.164  0.914    -1.000     0.268  0.251 -0.914      0.999    
##             p(__,2)1:H
## ply(c__,2)1           
## ply(c__,2)2           
## tretmntClng           
## tretmntHtng           
## pl(__,2)1:C           
## pl(__,2)2:C           
## pl(__,2)1:H           
## pl(__,2)2:H -0.914
#anova(CEWL_LMM_besta, type = "3", ddf = "Kenward-Roger")

CEWL_LMM_bestb <- lmerTest::lmer(data = temp_time_series,
                          CEWL_g_m2h ~ 
                             poly(calibrated_dorsal_temp, 2) * treatment +
                             (1|date/lizard_ID))
summary(CEWL_LMM_bestb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(calibrated_dorsal_temp, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 83028.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7628 -0.6345 -0.0833  0.5160  8.3769 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept)  9.996   3.162   
##  date           (Intercept) 21.953   4.685   
##  Residual                    1.554   1.246   
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                                    Estimate Std. Error
## (Intercept)                                           3.248      2.336
## poly(calibrated_dorsal_temp, 2)1                   1534.609     67.430
## poly(calibrated_dorsal_temp, 2)2                  -2711.185     75.106
## treatmentCooling                                      3.529      1.386
## treatmentHeating                                     14.341      1.468
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling -1350.765     67.457
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling  2852.577     75.134
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating -1365.961     67.469
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating  2666.790     75.141
##                                                          df t value Pr(>|t|)
## (Intercept)                                           5.353   1.390   0.2194
## poly(calibrated_dorsal_temp, 2)1                  25237.527  22.759  < 2e-16
## poly(calibrated_dorsal_temp, 2)2                  25236.097 -36.098  < 2e-16
## treatmentCooling                                     30.609   2.546   0.0162
## treatmentHeating                                     30.174   9.766 7.46e-11
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling 25237.523 -20.024  < 2e-16
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling 25236.096  37.967  < 2e-16
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating 25237.517 -20.246  < 2e-16
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating 25236.103  35.490  < 2e-16
##                                                      
## (Intercept)                                          
## poly(calibrated_dorsal_temp, 2)1                  ***
## poly(calibrated_dorsal_temp, 2)2                  ***
## treatmentCooling                                  *  
## treatmentHeating                                  ***
## poly(calibrated_dorsal_temp, 2)1:treatmentCooling ***
## poly(calibrated_dorsal_temp, 2)2:treatmentCooling ***
## poly(calibrated_dorsal_temp, 2)1:treatmentHeating ***
## poly(calibrated_dorsal_temp, 2)2:treatmentHeating ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(c__,2)1 -0.162                                                        
## ply(c__,2)2  0.163 -0.946                                                 
## tretmntClng -0.327  0.273    -0.276                                       
## tretmntHtng -0.315  0.257    -0.260     0.507                             
## pl(__,2)1:C  0.162 -1.000     0.946    -0.273 -0.257                      
## pl(__,2)2:C -0.163  0.946    -1.000     0.276  0.259 -0.946               
## pl(__,2)1:H  0.162 -0.999     0.946    -0.273 -0.257  0.999     -0.945    
## pl(__,2)2:H -0.163  0.946    -1.000     0.276  0.259 -0.945      0.999    
##             p(__,2)1:H
## ply(c__,2)1           
## ply(c__,2)2           
## tretmntClng           
## tretmntHtng           
## pl(__,2)1:C           
## pl(__,2)2:C           
## pl(__,2)1:H           
## pl(__,2)2:H -0.945
#anova(CEWL_LMM_bestb, type = "3", ddf = "Kenward-Roger")

# double check sample sizes
temp_time_series %>%
  group_by(lizard_ID, treatment) %>%
  summarise(n()) %>%
  group_by(treatment) %>%
  summarise(n())
## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.
## # A tibble: 3 × 2
##   treatment `n()`
##   <fct>     <int>
## 1 Control      11
## 2 Cooling      12
## 3 Heating      10

These are the best models for CEWL ~ temperature.

Export:

broom.mixed::tidy(CEWL_LMM_besta) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/CEWL_best_model_clotemp.csv")

broom.mixed::tidy(CEWL_LMM_bestb) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/CEWL_best_model_dorstemp.csv")

CEWL (raw) ~ time LMM

NOTE: Start times vary from 97-170 seconds after starting time series.

Make a polynomial model first:

rates_LMM_01 <- lme4::lmer(data = CEWL_time_series,
                          CEWL_g_m2h ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
drop1(rates_LMM_01)
## Single term deletions
## 
## Model:
## CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##                             npar    AIC
## <none>                            96906
## poly(time_min, 2):treatment    4 111987

Also make a linear model version:

rates_LMM_02 <- lme4::lmer(data = CEWL_time_series,
                          CEWL_g_m2h ~ 
                             time_min * treatment +
                             (1|date/lizard_ID))

anova(rates_LMM_02, rates_LMM_01)
## refitting model(s) with ML (instead of REML)
## Data: CEWL_time_series
## Models:
## rates_LMM_02: CEWL_g_m2h ~ time_min * treatment + (1 | date/lizard_ID)
## rates_LMM_01: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC logLik deviance  Chisq Df Pr(>Chisq)    
## rates_LMM_02    9 100913 100987 -50447   100895                         
## rates_LMM_01   12  96906  97005 -48441    96882 4012.3  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Polynomial is much better than linear.

Also compare to log version:

rates_LMM_03 <- lme4::lmer(data = CEWL_time_series,
                          CEWL_g_m2h ~ 
                             log(time_min) * treatment +
                             (1|date/lizard_ID))

anova(rates_LMM_02, rates_LMM_03)
## refitting model(s) with ML (instead of REML)
## Data: CEWL_time_series
## Models:
## rates_LMM_02: CEWL_g_m2h ~ time_min * treatment + (1 | date/lizard_ID)
## rates_LMM_03: CEWL_g_m2h ~ log(time_min) * treatment + (1 | date/lizard_ID)
##              npar    AIC    BIC logLik deviance  Chisq Df Pr(>Chisq)
## rates_LMM_02    9 100913 100987 -50447   100895                     
## rates_LMM_03    9  97121  97196 -48552    97103 3791.4  0
anova(rates_LMM_01, rates_LMM_03)
## refitting model(s) with ML (instead of REML)
## Data: CEWL_time_series
## Models:
## rates_LMM_03: CEWL_g_m2h ~ log(time_min) * treatment + (1 | date/lizard_ID)
## rates_LMM_01: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##              npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## rates_LMM_03    9 97121 97196 -48552    97103                         
## rates_LMM_01   12 96906 97005 -48441    96882 220.96  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The polynomial model is best (better than both log and lm options).

Check linear regression assumptions:

plot(rates_LMM_01)

vif(rates_LMM_01)
##                                 GVIF Df GVIF^(1/(2*Df))
## poly(time_min, 2)           9.050897  2        1.734494
## treatment                   1.000002  2        1.000000
## poly(time_min, 2):treatment 9.050904  4        1.317002

Again, there is a really weird spike, but otherwise there’s no fanning or pattern.

re-run with p-values:

rates_LMM_best <- lmerTest::lmer(data = CEWL_time_series,
                          CEWL_g_m2h ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
summary(rates_LMM_best)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_g_m2h ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##    Data: CEWL_time_series
## 
## REML criterion at convergence: 96853.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3325 -0.5700 -0.0438  0.5178  8.0025 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept)  9.976   3.158   
##  date           (Intercept) 17.817   4.221   
##  Residual                    1.754   1.324   
## Number of obs: 28405, groups:  lizard_ID:date, 37; date, 5
## 
## Fixed effects:
##                                      Estimate Std. Error        df t value
## (Intercept)                            14.957      2.099     5.189   7.126
## poly(time_min, 2)1                   -137.974      2.309 28362.024 -59.765
## poly(time_min, 2)2                     51.100      2.322 28362.016  22.008
## treatmentCooling                       -8.473      1.297    30.005  -6.533
## treatmentHeating                        3.416      1.280    30.061   2.669
## poly(time_min, 2)1:treatmentCooling   -88.851      3.376 28362.571 -26.319
## poly(time_min, 2)2:treatmentCooling    83.367      3.337 28362.083  24.980
## poly(time_min, 2)1:treatmentHeating   301.172      3.201 28362.093  94.092
## poly(time_min, 2)2:treatmentHeating  -108.483      3.201 28362.048 -33.893
##                                     Pr(>|t|)    
## (Intercept)                         0.000723 ***
## poly(time_min, 2)1                   < 2e-16 ***
## poly(time_min, 2)2                   < 2e-16 ***
## treatmentCooling                    3.17e-07 ***
## treatmentHeating                    0.012151 *  
## poly(time_min, 2)1:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating  < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1  0.000                                                    
## ply(tm_,2)2  0.000 -0.067                                             
## tretmntClng -0.305  0.001    0.000                                    
## tretmntHtng -0.317  0.001    0.000    0.499                           
## ply(_,2)1:C  0.000 -0.684    0.046    0.000  0.000                    
## ply(_,2)2:C  0.000  0.047   -0.696    0.000  0.000  0.008             
## ply(_,2)1:H  0.000 -0.721    0.049    0.000  0.000  0.493    -0.034   
## ply(_,2)2:H  0.000  0.049   -0.725    0.000  0.000 -0.033     0.505   
##             p(_,2)1:H
## ply(tm_,2)1          
## ply(tm_,2)2          
## tretmntClng          
## tretmntHtng          
## ply(_,2)1:C          
## ply(_,2)2:C          
## ply(_,2)1:H          
## ply(_,2)2:H -0.029
#anova(rates_LMM_best, type = "3", ddf = "Kenward-Roger")

# double check sample sizes
CEWL_time_series %>%
  group_by(lizard_ID, treatment) %>%
  summarise(n()) %>%
  group_by(treatment) %>%
  summarise(n())
## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.
## # A tibble: 3 × 2
##   treatment `n()`
##   <fct>     <int>
## 1 Control      12
## 2 Cooling      12
## 3 Heating      13

Export:

broom.mixed::tidy(rates_LMM_best) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/CEWL_best_model_time.csv")

Body Temperature ~ Time

linear models:

body_temp_LMM_01a <- lme4::lmer(data = temp_time_series,
                          calibrated_cloacal_temp ~ 
                             time_min * treatment +
                             (1|date/lizard_ID))

body_temp_LMM_01b <- lme4::lmer(data = temp_time_series,
                          calibrated_dorsal_temp ~ 
                             time_min * treatment +
                             (1|date/lizard_ID))

polynomial models:

body_temp_LMM_02a <- lme4::lmer(data = temp_time_series,
                          calibrated_cloacal_temp ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
anova(body_temp_LMM_02a, body_temp_LMM_01a)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## body_temp_LMM_01a: calibrated_cloacal_temp ~ time_min * treatment + (1 | date/lizard_ID)
## body_temp_LMM_02a: calibrated_cloacal_temp ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##                   npar   AIC   BIC logLik deviance Chisq Df Pr(>Chisq)    
## body_temp_LMM_01a    9 63638 63712 -31810    63620                        
## body_temp_LMM_02a   12 41291 41388 -20633    41267 22354  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
body_temp_LMM_02b <- lme4::lmer(data = temp_time_series,
                          calibrated_dorsal_temp ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
anova(body_temp_LMM_02b, body_temp_LMM_01b)
## refitting model(s) with ML (instead of REML)
## Data: temp_time_series
## Models:
## body_temp_LMM_01b: calibrated_dorsal_temp ~ time_min * treatment + (1 | date/lizard_ID)
## body_temp_LMM_02b: calibrated_dorsal_temp ~ poly(time_min, 2) * treatment + (1 | date/lizard_ID)
##                   npar   AIC   BIC logLik deviance  Chisq Df Pr(>Chisq)    
## body_temp_LMM_01b    9 77531 77604 -38756    77513                         
## body_temp_LMM_02b   12 71434 71532 -35705    71410 6102.6  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Again, the polynomial models are MUCH better.

body_temp_LMM_besta <- lmerTest::lmer(data = temp_time_series,
                          calibrated_cloacal_temp ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
summary(body_temp_LMM_besta)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: calibrated_cloacal_temp ~ poly(time_min, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 41254.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8989 -0.4649 -0.0011  0.3942  6.3225 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 3.0570   1.7484  
##  date           (Intercept) 0.1856   0.4308  
##  Residual                   0.2967   0.5447  
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                       Estimate Std. Error         df  t value
## (Intercept)                            27.4181     0.5621    16.5563   48.780
## poly(time_min, 2)1                     52.0361     0.9376 25212.0138   55.500
## poly(time_min, 2)2                      6.0643     0.9435 25212.0089    6.428
## treatmentCooling                       -4.0686     0.7317    26.4065   -5.560
## treatmentHeating                       -0.7262     0.7697    27.6621   -0.944
## poly(time_min, 2)1:treatmentCooling  -764.1781     1.3382 25212.2800 -571.054
## poly(time_min, 2)2:treatmentCooling   161.4913     1.3255 25212.0407  121.838
## poly(time_min, 2)1:treatmentHeating   716.1235     1.3557 25212.0589  528.233
## poly(time_min, 2)2:treatmentHeating   -63.6826     1.3531 25212.0326  -47.063
##                                     Pr(>|t|)    
## (Intercept)                          < 2e-16 ***
## poly(time_min, 2)1                   < 2e-16 ***
## poly(time_min, 2)2                  1.32e-10 ***
## treatmentCooling                    7.32e-06 ***
## treatmentHeating                       0.354    
## poly(time_min, 2)1:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating  < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001                                                    
## ply(tm_,2)2  0.000 -0.077                                             
## tretmntClng -0.676  0.000    0.000                                    
## tretmntHtng -0.647  0.000    0.000    0.490                           
## ply(_,2)1:C  0.000 -0.701    0.054    0.000  0.000                    
## ply(_,2)2:C  0.000  0.055   -0.712    0.000  0.000 -0.005             
## ply(_,2)1:H  0.000 -0.692    0.053    0.000  0.000  0.485    -0.038   
## ply(_,2)2:H  0.000  0.054   -0.697    0.000  0.000 -0.038     0.496   
##             p(_,2)1:H
## ply(tm_,2)1          
## ply(tm_,2)2          
## tretmntClng          
## tretmntHtng          
## ply(_,2)1:C          
## ply(_,2)2:C          
## ply(_,2)1:H          
## ply(_,2)2:H -0.021
write.csv(broom.mixed::tidy(body_temp_LMM_besta), 
          "./Results_Statistics/temp_time_best_model_clotemp.csv")

body_temp_LMM_bestb <- lmerTest::lmer(data = temp_time_series,
                          calibrated_dorsal_temp ~ 
                             poly(time_min, 2) * treatment +
                             (1|date/lizard_ID))
summary(body_temp_LMM_bestb)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: calibrated_dorsal_temp ~ poly(time_min, 2) * treatment + (1 |  
##     date/lizard_ID)
##    Data: temp_time_series
## 
## REML criterion at convergence: 71392.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8879 -0.3478 -0.0225  0.3300  6.0970 
## 
## Random effects:
##  Groups         Name        Variance Std.Dev.
##  lizard_ID:date (Intercept) 1.5384   1.2403  
##  date           (Intercept) 0.3563   0.5969  
##  Residual                   0.9812   0.9906  
## Number of obs: 25251, groups:  lizard_ID:date, 33; date, 5
## 
## Fixed effects:
##                                       Estimate Std. Error         df  t value
## (Intercept)                            27.7236     0.4606    11.4180   60.191
## poly(time_min, 2)1                     48.0115     1.7051 25212.0897   28.158
## poly(time_min, 2)2                      0.6860     1.7158 25212.0585    0.400
## treatmentCooling                       -2.8157     0.5208    26.3868   -5.406
## treatmentHeating                       -1.5969     0.5509    27.0887   -2.899
## poly(time_min, 2)1:treatmentCooling  -697.0769     2.4336 25213.7698 -286.444
## poly(time_min, 2)2:treatmentCooling   133.9832     2.4104 25212.2645   55.585
## poly(time_min, 2)1:treatmentHeating   484.8482     2.4654 25212.3649  196.660
## poly(time_min, 2)2:treatmentHeating   -42.9964     2.4608 25212.2024  -17.473
##                                     Pr(>|t|)    
## (Intercept)                         1.19e-15 ***
## poly(time_min, 2)1                   < 2e-16 ***
## poly(time_min, 2)2                   0.68930    
## treatmentCooling                    1.10e-05 ***
## treatmentHeating                     0.00734 ** 
## poly(time_min, 2)1:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)2:treatmentCooling  < 2e-16 ***
## poly(time_min, 2)1:treatmentHeating  < 2e-16 ***
## poly(time_min, 2)2:treatmentHeating  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) pl(_,2)1 pl(_,2)2 trtmnC trtmnH p(_,2)1:C p(_,2)2:C
## ply(tm_,2)1 -0.001                                                    
## ply(tm_,2)2  0.000 -0.077                                             
## tretmntClng -0.584  0.001    0.000                                    
## tretmntHtng -0.560  0.001    0.000    0.480                           
## ply(_,2)1:C  0.001 -0.701    0.054    0.000 -0.001                    
## ply(_,2)2:C  0.000  0.055   -0.712    0.000  0.000 -0.005             
## ply(_,2)1:H  0.001 -0.692    0.053   -0.001  0.000  0.485    -0.038   
## ply(_,2)2:H  0.000  0.054   -0.697    0.000  0.000 -0.038     0.496   
##             p(_,2)1:H
## ply(tm_,2)1          
## ply(tm_,2)2          
## tretmntClng          
## tretmntHtng          
## ply(_,2)1:C          
## ply(_,2)2:C          
## ply(_,2)1:H          
## ply(_,2)2:H -0.021
write.csv(broom.mixed::tidy(body_temp_LMM_bestb), 
          "./Results_Statistics/temp_time_best_model_dorstemp.csv")

Rate of Change

We want to know why each lizard had a different rate of change of CEWL during the 15-minute temperature treatment and measurement period.

First, calculate rate of change for each lizard:

change_LM <- lm(data = CEWL_time_series,
                          CEWL_g_m2h ~ 
                             time_min*lizard_ID)

Now get the trend/slope for each individual and make into a variable in df.

lizard_coeffs <- data.frame(emtrends(change_LM, "lizard_ID", var = "time_min")) %>%
   left_join(read_rds("Data/collection_dat_formatted.RDS"), 
             by = c("lizard_ID" = "individual_ID")) %>%
   dplyr::select(lizard_ID, rate_change = time_min.trend, 
                 SE, df, lower.CL, upper.CL,
                 date,
                 treatment,
                 mass_g, SVL_mm,
                 chamber_time_elapsed_hr
                 ) %>%
   dplyr::mutate(date_factor = as.factor(date))
summary(lizard_coeffs)
##    lizard_ID   rate_change            SE                 df       
##  401    : 1   Min.   :-0.8023   Min.   :0.008416   Min.   :28331  
##  402    : 1   1st Qu.:-0.3728   1st Qu.:0.008672   1st Qu.:28331  
##  404    : 1   Median :-0.1621   Median :0.008739   Median :28331  
##  406    : 1   Mean   :-0.1074   Mean   :0.009524   Mean   :28331  
##  407    : 1   3rd Qu.: 0.1088   3rd Qu.:0.008920   3rd Qu.:28331  
##  408    : 1   Max.   : 0.6960   Max.   :0.024080   Max.   :28331  
##  (Other):31                                                       
##     lower.CL           upper.CL             date              treatment 
##  Min.   :-0.81947   Min.   :-0.78520   Min.   :2021-10-05   Control:12  
##  1st Qu.:-0.38947   1st Qu.:-0.35612   1st Qu.:2021-10-11   Cooling:12  
##  Median :-0.17933   Median :-0.14491   Median :2021-10-12   Heating:13  
##  Mean   :-0.12610   Mean   :-0.08877   Mean   :2021-10-13               
##  3rd Qu.: 0.09151   3rd Qu.: 0.12617   3rd Qu.:2021-10-18               
##  Max.   : 0.67913   Max.   : 0.71289   Max.   :2021-10-19               
##                                                                         
##      mass_g          SVL_mm      chamber_time_elapsed_hr     date_factor
##  Min.   : 8.70   Min.   :60.00   Min.   :0.9333          2021-10-05:6   
##  1st Qu.:10.00   1st Qu.:64.00   1st Qu.:1.0000          2021-10-11:8   
##  Median :11.30   Median :66.00   Median :1.0333          2021-10-12:7   
##  Mean   :11.42   Mean   :66.24   Mean   :1.1005          2021-10-18:8   
##  3rd Qu.:12.60   3rd Qu.:69.00   3rd Qu.:1.1500          2021-10-19:8   
##  Max.   :15.30   Max.   :75.00   Max.   :1.5333                         
## 
length(unique(CEWL_time_series$lizard_ID))
## [1] 37

NOW we can look at how the RATE of CEWL change (CEWL change per MINUTE) is different among lizards.

change_LMM_1 <- lme4::lmer(data = lizard_coeffs,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                       mass_g + SVL_mm + 
                       chamber_time_elapsed_hr + 
                  # include date as random effect bc sig diff weather
                       (1|date_factor))

drop1(change_LMM_1)
## Single term deletions
## 
## Model:
## rate_change ~ treatment + mass_g + SVL_mm + chamber_time_elapsed_hr + 
##     (1 | date_factor)
##                         npar    AIC
## <none>                       21.929
## treatment                  2 42.116
## mass_g                     1 20.523
## SVL_mm                     1 21.118
## chamber_time_elapsed_hr    1 20.291
anova(change_LMM_1)
## Analysis of Variance Table
##                         npar  Sum Sq Mean Sq F value
## treatment                  2 2.56693 1.28346 17.8977
## mass_g                     1 0.00639 0.00639  0.0891
## SVL_mm                     1 0.08371 0.08371  1.1673
## chamber_time_elapsed_hr    1 0.02129 0.02129  0.2969

Drop mass first based on low SS and resulting AIC:

change_LMM_2 <- lme4::lmer(data = lizard_coeffs,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                       SVL_mm + 
                       chamber_time_elapsed_hr + 
                  # include date as random effect bc sig diff weather
                       (1|date_factor))
drop1(change_LMM_2)
## Single term deletions
## 
## Model:
## rate_change ~ treatment + SVL_mm + chamber_time_elapsed_hr + 
##     (1 | date_factor)
##                         npar    AIC
## <none>                       20.523
## treatment                  2 40.292
## SVL_mm                     1 19.118
## chamber_time_elapsed_hr    1 19.059
anova(change_LMM_2)
## Analysis of Variance Table
##                         npar  Sum Sq Mean Sq F value
## treatment                  2 2.56389 1.28195 17.5558
## SVL_mm                     1 0.02860 0.02860  0.3917
## chamber_time_elapsed_hr    1 0.03428 0.03428  0.4695

drop SVL:

change_LMM_3 <- lme4::lmer(data = lizard_coeffs,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                             chamber_time_elapsed_hr + 
                  # include date as random effect bc sig diff weather
                       (1|date_factor))
drop1(change_LMM_3)
## Single term deletions
## 
## Model:
## rate_change ~ treatment + chamber_time_elapsed_hr + (1 | date_factor)
##                         npar    AIC
## <none>                       19.118
## treatment                  2 38.317
## chamber_time_elapsed_hr    1 17.510
anova(change_LMM_3)
## Analysis of Variance Table
##                         npar  Sum Sq Mean Sq F value
## treatment                  2 2.56421 1.28211 17.8795
## chamber_time_elapsed_hr    1 0.02549 0.02549  0.3554

drop chamber_time_elapsed_hr:

change_LMM_4 <- lme4::lmer(data = lizard_coeffs,
                          rate_change ~ 
                             treatment + 
                       (1|date_factor))
drop1(change_LMM_4)
## boundary (singular) fit: see help('isSingular')
## Single term deletions
## 
## Model:
## rate_change ~ treatment + (1 | date_factor)
##           npar    AIC
## <none>         17.510
## treatment    2 38.876
anova(change_LMM_4)
## Analysis of Variance Table
##           npar Sum Sq Mean Sq F value
## treatment    2 2.5645  1.2823  18.286

The simplest model is our best/final model.

Check assumptions:

plot(change_LMM_4)

Assumptions look good.

Save with p-values:

change_LMM_best <- lmerTest::lmer(data = lizard_coeffs,
                          rate_change ~ # remember estimate is CEWL PER MINUTE
                             treatment + 
                       (1|date_factor))
summary(change_LMM_best)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rate_change ~ treatment + (1 | date_factor)
##    Data: lizard_coeffs
## 
## REML criterion at convergence: 16.9
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.94436 -0.62201 -0.04392  0.39606  1.77788 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  date_factor (Intercept) 0.01179  0.1086  
##  Residual                0.07012  0.2648  
## Number of obs: 37, groups:  date_factor, 5
## 
## Fixed effects:
##                  Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)       -0.2093     0.0908 12.1290  -2.305 0.039599 *  
## treatmentCooling  -0.1793     0.1085 29.9201  -1.653 0.108748    
## treatmentHeating   0.4475     0.1068 30.3121   4.191 0.000222 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trtmnC
## tretmntClng -0.593       
## tretmntHtng -0.611  0.504
anova(change_LMM_best, type = "3", ddf = "Kenward-Roger")
## Type III Analysis of Variance Table with Kenward-Roger's method
##           Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## treatment 2.5102  1.2551     2 30.708  17.898 7.034e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Treatment group explained a significant amount of the variation (F(3,30)=18.29, p < 0.0001).

Export:

broom.mixed::tidy(change_LMM_best) %>%
  mutate(estimate = round(estimate, digits = 2),
         std.error = round(std.error, digits = 2),
         statistic = round(statistic, digits = 2),
         df = round(df, digits = 0)) %>%
  dplyr::select(effect, group, term, estimate, 
                std.error, statistic, df, p.value) %>%
  write.csv("./Results_Statistics/change_CEWL_best_model.csv")

also get pairwise differences for change ~ tmt, which I will use to make the boxplot later:

emmeans <- data.frame(emmeans(change_LMM_best, "treatment"))
emmeans
##   treatment     emmean         SE       df    lower.CL    upper.CL
## 1   Control -0.2093260 0.09102635 13.09716 -0.40582833 -0.01282375
## 2   Cooling -0.3886322 0.09189944 12.94683 -0.58725181 -0.19001265
## 3   Heating  0.2381211 0.08845166 11.97139  0.04535039  0.43089180
emmean_ps <- data.frame(pairs(emmeans(change_LMM_best, "treatment"))) # p values
emmean_ps
##            contrast   estimate        SE       df   t.ratio      p.value
## 1 Control - Cooling  0.1793062 0.1087886 30.33842  1.648208 2.414730e-01
## 2 Control - Heating -0.4474471 0.1074827 30.69471 -4.162968 6.708347e-04
## 3 Cooling - Heating -0.6267533 0.1082740 31.09322 -5.788584 6.511073e-06

Figures

COLOR

cool <- c(brewer.pal(11, "RdBu")[c(11)])
heat <- c(brewer.pal(11, "RdBu")[c(2)])
control <- c(brewer.pal(9, "Greys")[c(5)])
cntrls <- c(brewer.pal(9, "Greys")[c(7, 5, 3, 5, 4, 5, 3, 6, 7, 3, 6)])

cloacal temp ~ time

ggplot(temp_time_series) +
   aes(x = time_sec, 
       y = calibrated_cloacal_temp,
       color = treatment) +
   geom_point(size = 0.01, 
              alpha = 0.2) +
   geom_smooth(method = "loess", 
               se = F, 
               size = 2,
               na.rm = TRUE) +
   theme_classic() +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time (seconds)" , 
        y = "Internal Body Temperature (°C)", 
        color = "Treatment", se = FALSE) -> ctemp_time
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
ctemp_time
## `geom_smooth()` using formula = 'y ~ x'

dorsal temp ~ time

ggplot(temp_time_series) +
   aes(x = time_sec, 
       y = (calibrated_dorsal_temp),
       color = treatment) +
   geom_point(size = 0.01, 
              alpha = 0.2) +
   geom_smooth(method = "loess", 
               se = F, 
               size = 2,
               na.rm = TRUE) +
   theme_classic() +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x= "Time (seconds)" , 
        y= "Surface Body Temperature (°C)", 
        color ="Treatment") -> dtemp_time
dtemp_time
## `geom_smooth()` using formula = 'y ~ x'

CEWL (raw) ~ time

ggplot(CEWL_time_series) +
   aes(x= time_sec, 
       y = CEWL_g_m2h,
       color = treatment) +
   geom_point(size = 0.01, 
              alpha = 0.25) +
   geom_smooth(method = "loess", 
               size = 2.5, 
               se = F) +
   theme_classic() +
   theme(text = element_text(size = 15)) + 
   theme(axis.text.x = element_text(size = 10)) +
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time (seconds)" , 
        y = bquote('CEWL (g/'*m^2*'h)'), 
        color = "Treatment", 
        se = FALSE) +
   scale_x_continuous(limits = c(0, 925), 
                      breaks=seq(0, 950, 100)) -> CEWL_time
CEWL_time
## `geom_smooth()` using formula = 'y ~ x'

legend

This code uses one of our plots to make a large, detailed legend, which we can save and use in ggarrange later.

ggplot(temp_time_series) +
   
   geom_smooth(aes(x = (calibrated_cloacal_temp), 
       y = CEWL_g_m2h,
       color = treatment),
       method = "lm", 
        formula = y ~ poly(x, 2),
               size = 1, 
               se = F) +
  geom_point(aes(x = (calibrated_cloacal_temp), 
       y = CEWL_g_m2h,
       fill = treatment,
       shape = treatment),
       size = 4, 
       alpha = 1) +
  
  scale_x_continuous(limits = c(14, 39),
                     breaks = c(seq(15, 35, 5)),
                     labels = c(seq(15, 35, 5))) +
  theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = c(0.5,0.5)) +
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") +
   scale_fill_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  scale_color_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment")  -> fancy_legend
fancy_legend # view

# save
LEGEND <- as_ggplot(get_legend(fancy_legend)) 

# export
ggsave(filename = "legend_only.pdf",
       plot = LEGEND,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 40, height = 30)

CEWL (raw) ~ cloacal temperature

ggplot(temp_time_series) +
   geom_point(aes(x = (calibrated_cloacal_temp), 
       y = CEWL_g_m2h,
       shape = treatment,
       color = treatment),
       size = 0.01, 
               alpha = 0.1) +
   geom_smooth(aes(x = (calibrated_cloacal_temp), 
       y = CEWL_g_m2h,
       color = treatment),
       method = "lm", 
        formula = y ~ poly(x, 2),
               size = 1.5, 
               se = F) +
  scale_x_continuous(limits = c(12.5, 40),
                     breaks = c(seq(15, 35, 5)),
                     labels = c(seq(15, 35, 5))) +

  # annotations for cooling group
   annotate(geom = "point", x = 14.3, y = 5, 
           size = 4, pch = 25, fill = cool) + 
   annotate("text", x = 12.6, y = 5, label = "T[15]", parse = T,
           size = 5, family = "sans", color = cool) +
   annotate(geom = "point", x = 35.2, y = 11.8, 
           size = 4, pch = 25, fill = cool) + 
   annotate("text", x = 36.7, y = 11.8, label = "T[1]", parse = T,
           size = 5, family = "sans", color = cool) +
  
  # annotations for heating group
   annotate(geom = "point", x = 15.9, y = 14, 
           size = 4, pch = 24, fill = heat) +  
   annotate("text", x = 14.3, y = 14, label = "T[1]", parse = T,
           size = 5, family = "sans", color = heat) +
   annotate(geom = "point", x = 38.3, y = 25.2, 
           size = 4, pch = 24, fill = heat) +
   annotate("text", x = 40, y = 25.2, label = "T[15]", parse = T,
           size = 5, family = "sans", color = heat) +
   
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
  
   scale_color_manual(values = c("Control" = control,
                                   "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
  
   labs(x = "Internal Body Temperature (°C)" , 
        y = bquote('CEWL (g '*m^-2*' '*h^-1*')'), 
        se = FALSE) -> CEWL_ctemp
CEWL_ctemp

# save
ggsave(filename = "CEWL_int_body_temp.pdf",
       plot = CEWL_ctemp,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 150, height = 100)

CEWL Control

# get control data only
temp_time_series %>%
  dplyr::filter(treatment == "Control") %>%
  # pipe into ggplot
  ggplot() +
   aes(x = calibrated_cloacal_temp, 
       y = CEWL_g_m2h,
       color = lizard_ID) +
  
   geom_point(size = 0.1, 
               alpha = 0.1) +
   geom_smooth(method = "lm", 
               size = 1, 
               se = F) +

   ylim(0,32) +
   xlim(25,31) +
  
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.text.align = 0,
        legend.position = "none",
        plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
  
   scale_color_manual(values = cntrls) +
  
   labs(x = "Internal Body Temperature (°C)" , 
        y = bquote('CEWL (g '*m^-2*' '*h^-1*')'), 
        color = "Treatment", 
        se = FALSE) -> CEWL_control
CEWL_control
## `geom_smooth()` using formula = 'y ~ x'

# save
ggsave(filename = "CEWL_temp_control_lm.pdf",
       plot = CEWL_control,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)
## `geom_smooth()` using formula = 'y ~ x'

CEWL (relative) ~ time

# check n lizard sample size
CEWL_time_series %>%
  dplyr::filter(complete.cases(time_min, relative_CEWL)) %>%
  group_by(lizard_ID) %>%
  summarise(max(relative_CEWL)) %>%
  nrow()
## [1] 37
# plot
ggplot() +
   geom_point(data = CEWL_time_series,
              aes(x= time_min, 
                  y = relative_CEWL,
                  shape = treatment,
                  color = treatment),
              size = 0.01, 
              alpha = 0.1) +
  
   geom_hline(yintercept = 0, size = 0.5,
              linetype = "dashed", color = "black") +
   
   geom_smooth(data = CEWL_time_series,
               aes(x = time_min, 
                   y = relative_CEWL,
                   color = treatment),
               method = "lm", 
               formula = y ~ poly((x), 2),
               size = 1.5, 
               se = F) +
   
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none") +
  
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat),
                      name = "Treatment") +
  
  annotate(geom = "point", x = 15, y = -3.9, 
           size = 4, pch = 21, fill = control) +
   annotate(geom = "point", x = 15, y = -7.5, 
           size = 4, pch = 25, fill = cool) + 
   annotate(geom = "point", x = 15, y = 3, 
           size = 4, pch = 24, fill = heat) +
  scale_x_continuous(limits = c(1, 15), 
                      breaks = seq(1, 15, 2)) +

   labs(x = "Time (minutes)" , 
        y = bquote('CEWL (g '*m^-2*' '*h^-1*')')) -> CEWLr_time
  
CEWLr_time
## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 4 rows containing missing values (`geom_point()`).

# save
ggsave(filename = "CEWL_relative_time.pdf",
       plot = CEWLr_time,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 150, height = 100)
## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).
## Removed 4 rows containing missing values (`geom_point()`).

CEWL Rate of Change

ggplot() +
  geom_jitter(data = lizard_coeffs,
              aes(x = treatment,
                   y = rate_change,
                   shape = treatment,
                   color = treatment,
                   fill = treatment),
                   size = 0.8,
               alpha = 0.3,
              position = position_jitter(height = 0, width = 0.2)) +
  geom_errorbar(data = emmeans,
                aes(x = treatment,
                    y = emmean, 
                    color = treatment,
                    group = treatment,
                    ymin = lower.CL, 
                    ymax = upper.CL),
                width = .1,
                alpha = 0.9) +
  geom_point(data = emmeans, 
             aes(x = treatment,
                   y = emmean, 
                 shape = treatment,
                 fill = treatment), 
              color = "black",
              size = 4) +
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none") +
  
   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
  scale_fill_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
  
  scale_shape_manual(values = c(21,25,24),
                     name = "Treatment") + 
  
   geom_hline(yintercept = 0, size = 0.5,
              linetype = "dashed", color = "black") +
  
   annotate("text", x = 1, y = 0.42, label = "A", size = 3) +
   annotate("text", x = 2, y = 0.3, label = "A", size = 3) +
   annotate("text", x = 3, y = 0.9, label = "B", size = 3) +
  
   labs(x = "" , 
        y = expression(Delta ~ 'CEWL '*min^-1*''),
        color = "Treatment") -> rate_change_CEWL_fig
rate_change_CEWL_fig

# save
ggsave(filename = "CEWL_change_min.pdf",
       plot = rate_change_CEWL_fig,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)

Experimental Design Skematic

We created reference values to use to make a figure that would show each step of our experiment.

exp_design <- read.csv("./Data/exp design.csv", 
                            header = TRUE,
                            fileEncoding="UTF-8-BOM")
ggplot(exp_design) +
   aes(x = time_min, 
       y = temp_C,
       color = treatment) +
     geom_line(size = 1) +
   geom_vline(xintercept = 60,
              alpha = 0.5,
              size = 0.5,
              linetype = "dashed", color = "black") +
   geom_text(aes(30, 37, family = "sans"), label = "a", color = "black", size = 3) +
   geom_vline(xintercept = 65,
              alpha = 0.5,
              size = 0.5,
              linetype = "dashed", color = "black") +
   geom_text(aes(62.5, 37, family = "sans"), label = "b", color = "black", size = 3) +
   geom_vline(xintercept = 80,
              alpha = 0.5,
              size = 0.5,
              linetype = "dashed", color = "black") +
   geom_text(aes(72.5, 37, family = "sans"), label = "c", color = "black", size = 3) +
   theme_classic() +
   theme(text = element_text(color = "black", 
                            family = "sans", 
                            size = 12),
        axis.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 8),
        legend.text = element_text(color = "black", 
                                 family = "sans", 
                                 size = 12),
        legend.position = "none",
        axis.text.x = element_blank(),
        #axis.line.x = element_blank(),
        axis.ticks.x = element_blank()) +

   scale_color_manual(values = c("Control" = control,
                                     "Cooling" = cool,  
                                     "Heating" = heat)) +
   labs(x = "Time" , 
        y = "Body Temperature (°C)", 
        color = "Treatment") +
   
    annotate(geom = "point", x = 80, y = 25, 
           size = 3, pch = 21, fill = control) +
   annotate(geom = "point", x = 80, y = 20, 
           size = 3, pch = 25, fill = cool) + 
   annotate(geom = "point", x = 80, y = 30, 
           size = 3, pch = 24, fill = heat) +
   
   scale_y_continuous(limits = c(14, 37), 
                      breaks=seq(15, 35, 5)) -> methods_schmatic_fig
methods_schmatic_fig

# save
ggsave(filename = "methods_schematic.pdf",
       plot = methods_schmatic_fig,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 80, height = 60)

Save Grouped Figures

Relative CEWL & Rate of Change

multi_fig_relative <- ggarrange(CEWLr_time, 
          ggarrange(rate_change_CEWL_fig, LEGEND, 
                    nrow = 1, ncol = 2,
                    widths = c(2,1),
                    align = "hv"
                    ),
          nrow = 2, ncol = 1,
          heights = c(2,1),
          labels = c("A", "B"))
## Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 4 rows containing missing values (`geom_point()`).
multi_fig_relative

ggsave(filename = "CEWL_time_multi_fig.pdf",
       plot = multi_fig_relative,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 120, height = 160)

Raw CEWL ~ Cloacal temp & control zoom-in

multi_fig_CEWL_by_temp <- ggarrange(CEWL_ctemp, 
                  ggarrange(CEWL_control, LEGEND, 
                            nrow = 1, ncol = 2,
                            widths = c(2, 1),
                            align = "hv"
                            ),
                  nrow = 2, ncol = 1,
                  heights = c(2,1),
          labels = c("A", "B"))
## `geom_smooth()` using formula = 'y ~ x'
multi_fig_CEWL_by_temp

ggsave(filename = "CEWL_temp_multi_fig.pdf",
       plot = multi_fig_CEWL_by_temp,
       path = "./Results_Figures",
       device = "pdf",
       dpi = 600,
       units = "mm",
       width = 120, height = 160)